Metalenses phase characterization by multi-distance phase retrieval

Metalens, characterized by their unique functions and distinctive physical properties, have gained significant attention for their potential applications. To further optimize the performance of metalens, it is necessary to characterize the phase modulation of the metalens. In this study, we present a multi-distance phase retrieval system based on optical field scanning and discuss its convergence and robustness. Our findings indicate that the system is capable of retrieving the phase distribution of the metalens as long as the measurement noise is low and the total length of the scanned light field is sufficiently long. This system enables the analysis of focal length and aberration by utilizing the computed phase distribution. We extend our investigation to measure the phase distribution of the metalens operating in the near-infrared (NIR) spectrum and identify the impact of defects in the sample on the phase. Additionally, we conduct a comparative analysis of the phase distribution of the metalens in air and ethanol and observe the variations in the phase modulation of the metalens in different working mediums. Our system provides a straightforward method for the phase characterization of metalens, aiding in optimizing the metalens design and functionality.


Phase retrieval algorithm for MDPR system
Initially, generate a random sample-plane field distribution  , , where  is the sample plane,  is the number of iterations.Use the angular spectrum propagation method to calculate the field distribution corresponding to the given sample-plane field distribution at that position   .
Extract the phase of the calculated field distribution    at position   , combine it with the measured intensity distribution   , and obtain the updated far-field field distribution.To enhance the iterative convergence of the algorithm, double-feedback introduces data from the k-1, k-2 iterations when  > 2  ,+1 = (1 +  + ) ′ , −  ,−1 −  ,−2 (6) Through experimentation with actual data, it was found that parameters a and b exhibit the best convergence when their values are between 0.7 and 0.75.
The 'Double-feedback' step (6) in the algorithm is primarily aimed at enhancing the    In practical industrial applications, the convergence speed of the iterative algorithm is crucial.Therefore, we need to introduce more information to assist the algorithm in converging faster.Compared to traditional lensless imaging systems, our scanning system, based on a microscopic imaging system, allows for the capture of intensity distribution at the sample plane.Inspired by Kalman filtering, we employed a method of controlling weights ω to achieve balancing the intensity distribution obtained from measurements with the intensity distribution obtained through iteration in Eq.S5 = ( ′′ )) ( 7) By adjusting weights ω, we balance the direct measurement results (intensity distribution at the sample plane) with the indirect measurement results (intensity distribution obtained during the algorithm iteration process).This introduces measurement information at the sample plane, making it closer to the target value and allowing the algorithm to converge faster.Qualitatively, as the iteration progresses, we know that the iterative result will continuously approach the true value.Therefore, we choose to balance the iterative field intensity at the sample plane with the measured intensity using linearly changing weights ω from 1 to 0.

Fig.S3
shows the variation of MSE to the number of iterations across varying scan lengths, depicting both the algorithm's performance in its raw form and its enhancement with Kalman filtering, as depicted in Fig. S3(a) and Fig. S3(b), respectively.After introducing the Kalman filtering step, the iteration was accelerated by 1/3, reducing from 15 to 10 iterations.Additionally, with the introduction of sample intensity, the MSE decreases monotonically during the iteration process, indicating that the iterative algorithm is more stable.Due to the focusing effect of the metalens on the light beam, the intensity distribution behind the sample undergoes drastic changes.One case is when the scanning length is compatible with the focal length.When we require the optical path to capture clear images near the sample surface, images closer to the focus may be overexposed.
Another case is that under immersion conditions, multiple interfaces of liquid, container walls, and air reflections introduce coherent noise into the optical path, typically appearing as bright speckles.The intensity measurement here is erroneous, so such erroneous intensity information should not be introduced into the iteration process.For these reasons, we introduce a mask function M in updating the far-field intensity in Eq.S3 and update it as Where is the Mask function for the selected region.By introducing  , we address the impact of coherent noise in the optical path under immersion conditions, as shown in All the optimizations we discussed so far are based on optimizing the noise of the measurement system itself, which can accelerate the convergence speed of the algorithm and enhance its robustness.However, in actual experimental processes, we found that occasional errors during measurements can greatly affect the convergence of the algorithm.These occasional errors may arise from the optical platform's vibration or the introduction of ambient light.Especially in immersion conditions, vibrations from other instruments on the same optical platform can affect the smoothness of the liquid surface.External occasional disturbances can lead to significant errors in the intensity or position of an image during the scanning process, affecting the convergence and accuracy of the algorithm.Hence, it's imperative to identify images affected by these incidental disturbances and exclude them from the phase retrieval process.
We realized that during the iterations, we achieved a balance across intensity measurement in all positions at the sample plane.Since incidental disturbances occur only in a few images, the corresponding sample surface field distributions for these images should exhibit significant differences compared to those for images unaffected by the disturbances.By calculating the field distributions between the sample surfaces corresponding to different images, we can eliminate the images affected by disturbances.
After the iterations have stabilized, compute the RM deviation between    .
A simulated result of the RMS deviation distribution is shown in Fig. S6.We introduced arbitrary disturbances to the measurements of position and intensity at five locations in the simulation, causing them to deviate from the overall measurement results.Our method can filter out images corresponding to larger RMS values numbered 2, 9, 22, 23, and 27, which also means that after removing these images, the system can retrieve the phase precisely.The time complexity of the phase retrieval algorithm needs to be considered in two parts: the time complexity of the iteration process and the time complexity of angular spectrum propagation.For iterating T times using  images, the time complexity is M×T.The time complexity of the angular spectrum propagation process mainly comes from the calculation of the k-space distribution of the optical field sized N using 2-D FFT, which has a time complexity of  2 log .Therefore, the total time complexity of the algorithm is M×T× 2 log .On a PC configured with AMD Ryzen 7 5800H with Radeon Graphics and 16GB RAM, the time for phase retrieval of the NIR sample using 30 images sized 256×256 for 100 iterations in non-parallel conditions is 13 seconds.
2. The accuracy of the MDPR system and its consistency with off-axis interferometry.

consistency with off-axis interferometry.
The consistency of measurement results across different methods is a key issue in metrology.We used off-axis interferometry and the MDPR system to measure the phase of the same metalens.The difference in aberration between the two methods is shown in Fig. S7(a), with a RMS deviation of 0.16.

The accuracy of the MDPR system
We calibrated the phase accuracy of our method using standard samples, and the result is shown in Fig. S8.We scanned a pinhole in 532 nm without internal phase modulation inside, so all phase fluctuations in the measurement result are due to measurement errors.Using a 3 mm scan field data, the standard deviation of the phase retrieval obtained is 0.02 rad, which corresponds to a path length error of 1.6 nm (0.3% λ).

The scanning field measurement optical path for metalens in different mediums
The beam size of the incident light beam is expanded by a beam expander.The beam size is about 3 mm, which is sufficiently large compared to the sample size.Therefore, it can be assumed that the incident light in the sample region is collimated.Thus, the measured phase distribution is entirely due to the modulation of the sample.
the metalens is placed on the wall of the transparent water tank.The light field scanning measurement system is shown in Fig. S16.We placed the metalens in the air (Fig. S16 (a)) and immersed it in ethanol absolute (Fig. S16 (b)), respectively.The light field intensity distribution of the metalens was collected in these two environments.

The influence of scan field length on phase reproduction results
The mean square error does not monotonically decrease with the increase of , But there is a threshold instead.This is because before reaching the threshold, the phase distribution obtained by the algorithm differs significantly from the standard value, and at this point, the MSE cannot be used to measure the difference between them quantitatively.
The convergence behavior of the algorithm concerning  can be understood as follows: due to the increasing angle of deflection of light at different positions in a metalens with increasing distance from the lens center, the spatial frequency modulation of the light field by the metalens also increases with distance from the lens center.As shown in

Phase retrieval for other metasurfaces
Since the implementation of all-optical diffractive neural networks in the THZ regime (Science 361(6406), 1004(Science 361(6406), -1008(Science 361(6406), , 2018)), optical neural networks have garnered widespread attention due to their high computational speed, low power consumption, and minimal heat generation.To achieve device miniaturization, there is a collective aspiration to imply optical neural network architectures by metasurfaces within the visible light spectrum.As optical neural network chips shrink in size (on the order of hundreds of micrometers), alignment issues for the neural network chip emerge.
Deviations in the positions between layers can reduce the accuracy of neural networks and increase system energy consumption, thus limiting the increase in the number of layers in optical neural networks operating in the visible light spectrum, while increasing the number of layers remains advantageous for optical neural networks to enhance accuracy and realize complex functionalities (Light: Science & Applications, 10 ), 196, 2021).At the same time, for each layer of metasurface optical neural network chips, there will inevitably be processing errors.The processing and alignment errors will couple, and it is impossible to directly calibrate the alignment errors of multilayer optical neural networks from a functional perspective.It is necessary to measure the overall optical modulation of the multilayer network and the phase modulation of each layer to decouple the alignment errors of the system from the processing errors.For this section, a two-layer metasurface-based optical neural network is used as an example to present the application scenarios of this phase measurement technology.As  Here, we demonstrate a simulation-based alignment process of a double-layer metasurfaces system.During the alignment process, we first align the layer-1 with the optical path.By scanning the intensity distribution, the position of the layer-1 can be obtained, and the layer-1 can be adjusted to the center of the optical path, as shown in updated N far-field light field distributions are backward propagated to the sample plane position.   =  −1  −    ( ′   ) (4)To balance the information across the various intensity distributions, the obtained N sample-plane field distributions are averaged.
algorithm's noise tolerance and convergence speed.The selection of parameters a and b is influenced by the measurement conditions, especially influenced by the noise level and translation stage accuracy during measurements.We simulated the distribution of MSE between the phase retrieval results and the true values under different values of a and b for the same measurement conditions, as shown in Fig.S1.It can be observed that within a certain range (depicted by the white region in Fig.S1), the values corresponding to a and b can lead the algorithm to converge near the true values.

Fig. S1 :
Fig. S1: The sampled results of the same metalens phase distribution at different sampling intervals.The noise level is 0.1, and the uncertainty of the displacement stage is 1 µm after 80 iterations.

Fig. S2
Fig.S2illustrates the relationship between the MSE and the total scanning field length and iteration times.It can be observed that when the scanning field length is insufficient, the MSE remains unchanged with iteration number increases, indicating the algorithm is trapped in a local minimum.However, when the scanning field length is sufficiently long, the algorithm converges to the true value within ten iterations, avoiding the local minimum.

Fig. S2 :
Fig. S2: Relationship between RMS error, iteration number, and scanning field length with 30 pictures used for phase retrieval.

Fig. S3 :
Fig. S3: The variation of MSE with the number of iterations and the total length of the scan.(a) Without Kalman filtering (b) With Kalman filtering.

Fig. S6 :
Fig.S6: The near-field field distributions corresponding to different farfield positions.,  are the image numbers.

Fig. S8 :
Fig. S8: Phase retrieval results for the standard sample. is the normalization coefficient for length.

Fig. S9 :
Fig. S9: The measured spectrum at a laser wavelength of 1560 nm.

Fig. S12
Fig.S12 shows the sample image during the scanning.Due to the sample relying on cross-polarization for field focusing, non-sample regions exhibit dark noise after polarization separation, which is zero after subtracting the background.

Fig
Fig. S11.The optical image (A) and the SEM image (B) of the fabricated NIR metalens.

Fig. S14 :
Fig. S14: Statistical data of the noise

Fig. S15 :
Fig. S15: Zernike expansion coefficients for wave aberrations from different sources.(a) Distribution of wave aberrations due to processing and design errors.(b) Zernike expansion of aberration in (a).(c) Total aberration of the sample.(d)-(f) Zernike expansion of aberrations from hole(region ③),dust(region ①),and stain(region ④) respectively.

Fig. S16 :
Fig. S16: The light field scanning measurement system.(a) The metalens is in the air.(b) The metalens is immersed in ethanol absolute.

Fig. S17 :
Fig. S17: The effect of changing the working medium on the sample.(a) The modulation of the meta-atom phase in different media as a function of radius (b) The discretized phase distribution profile of metalens.The black line is the phase distribution in air.The red line is the phase distribution in alcohol.

Fig
Fig.S19.(c), for regions with high-frequency optical modulation at the edge of the metalens, the corresponding light field distribution changes dramatically in the zdirection.Rapidly changing intensity in z-direction implies more information, requiring only a short scanning distance to obtain sufficient information.For regions with lowfrequency optical modulation at the center of the metalens, the corresponding light field distribution changes smoothly in the z-direction, implying less information and requiring a longer scanning distance to obtain sufficient information.After the scanning total length reaches a threshold, increasing the scanning length only provides redundant information, which helps the robustness of the algorithm but does not affect the convergence.Fig.S19.(d) and (e) show the retrieval results before and after the threshold, indicates before the threshold only regions with high frequency can converge.For the number of images increasing the number of images also introduces more information.Considering both factors, the threshold forms a shape similar to a hyperbola in Fig.S19.(a) and (b).

Fig. S19 :
Fig. S19: MSE distribution of MDPR system in metalens measurement with sampling picture number N and distance between pictures L for sample working at 1550 nm with focal length (a) 0.6 mm and (b) 1 mm.(c) Schematic diagram of light field distribution corresponding to different frequency regions of the sample.

Figure S20 :
Figure S20: Phase retrieval results of NIR samples with different picture numbers.(a) Results with 4 pictures.(b)Results with 5 pictures.(c) Results with 10 pictures.

Figure S21 :
Figure S21: Alignment process diagram for bilayer optical neural network (ONN).(a) Flowchart for alignment.(b) Schematic diagram of ONN working principle and position offset.(c) Phase retrieval for layer-1.(d) Phase retrieval for bilayer ONN.(e) ONN after correcting lateral offset.

Fig. S21 .
Fig.S21.(b)shows,there is a certain displacement between two phase-modulating metasurfaces and the input and output layers due to alignment accuracy issues.Taking the layer-2 as an example, it could experience a displacement  relative to the optical axis in the x-y plane, along with a rotation angle .Simultaneously, there could be a deviation ∆ from the standard value in the z-direction relative to the spacing with the layer-1, while a 2% error in the z-direction position results in a 1% decrease in

Figure S22 :
Figure S22: Phase retrieval results for metasurfaces ONN chips.(a) The layer-1, (b) The layer-2 and its detail enlargement.(c) Phase distribution at the dashed line in (b).

Fig. S23 .
Fig.S23.(a).Afterward, we added the layer-2.At this point, the phase distribution is contributed by two parts: the diffraction from the layer-1 and the modulation from the layer-2.Fig.S23.(b) shows the phase modulation of the bilayer optical neural network system measured by MDPR.Therefore, it is difficult to directly determine the precise boundary position of the layer-2.Since the distance between layer-1 and layer-2 is much larger than the size of the metasurfaces (about 100 times), the contribution of diffraction from layer-1 at the position of layer-2 is low-frequency.By removing the contribution

Figure S23 :
Figure S23: Simulation results of the 2 layers optical neural network.(a) Place the layer-1 in the center of the optical path.(b)Phase retrieval results for the bilayer sample.(c)High-pass filtering results of (b).

Figure S24 :
Figure S24: Simulation results of the 2 layers optical neural network.(a) The phase distribution of the bilayer network after xy-plane alignment.(b) Remove the contribution of the layer-2 from (a).(c) the low-pass filtering of (b).(d) The standard deviation between (c) and the standard distribution.The lowest point position corresponds to the deviation of the obtained layer-2 position relative to the ideal position.